Efficient Computation of Dendritic Microstructures using Adaptive Mesh Refinement 
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We study dendritic microstructure evolution using an adaptive grid, finite element method applied 
to a phase-field model. The computational complexity of our algorithm, per unit time, scales 
linearly with system size, rather than the quadratic variation given by standard uniform mesh 
schemes. Time-dependent calculations in two dimensions are in good agreement with the predictions 
of solvability theory, and can be extended to three dimensions and small undercoolings. 
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Dendrites are the primary component of solidification 
microstructures in metals. The formation, shape, speed 
and size of dendritic microstructures has been a topic of 
intense study in the past 10-15 years. Experiments [QJ^] 
by Glicksman and coworkers on succinonitrile (SCN) and 
other transparent analogues of metals have been accurate 
enough to provide tests of theories of dendritic growth, 
and have stimulated considerable theoretical progress 
The experiments have clearly demonstrated that 
naturally growing dendrites possess a unique steady state 
tip, characterized by its velocity, radius of curvature and 
shape, which leads the to time-dependent sidebranched 
dendrite as it propagates. 

The earliest theories of dendritic growth solved for the 
diffusion field around a self-similar body of revolution 
propagating at constant speed ||;@]. I n these studies 
the diffusion field determines the product of the dendrite 
velocity and tip radius, but neither quantity by itself. 
Adding capillarity effects to the theory predicts a unique 
maximum growth speed, {| but experiments showed that 
this point does not represent the operating state for real 
dendrites. 

The goal of contemporary research has been to predict 
steady state features of dendritic growth and to com- 
pute time-dependent microstructures from numerical so- 
lutions of the equations of motion. The purpose of this 
letter is to present a computationally efficient method 
for time-dependent calculations, and to verify that the 
steady state properties are in excellent agreement with 
those predicted by analysis of the steady state problem. 

Insight into the steady state dendrite problem was first 
obtained from local models Q describing the evo- 

lution of the interface, and incorporating the features of 
the bulk phases into the governing equation of motion 
for the interface. These models were the first [|l0| to 
show that a nonzero dendrite velocity is obtained only 
if a source of anisotropy - for example, anisotropic in- 
terface energy - is present in the description of dendritic 
evolution. Subsequently, it was shown that the spectrum 
of allowed steady state velocities is discrete, rather than 
continuous, and the role of anisotropy was understood 
theoretically, both in the local models and the full mov- 



ing boundary problem |^,^,|lj,^9| . Moreover, only the 
fastest of a spectrum of steady state velocities is sta- 
ble, thus forming the operating state of the dendrite. 
It is widely believed that sidebranching is generated by 
thermal or other statistical fluctuations on a microscopic 
scale, which are amplified by advective diffusion. This 
body of theoretical work is generally known as solvabil- 
ity theory. 

Numerically solving the timc-dcpcndcnt Stefan prob- 
lem, or variations of it, is difficult, requiring front track- 
ing and lattice deformation to contain the interface at 
predefined locations on the grid [pcfl . These difficulties 
have been partially addressed by the introduction of the 
phase-field model, which introduces an auxiliary continu- 
ous order parameter <j)(r) that couples to the evolution of 
the thermal field. The phase field interpolates between 
the solid and liquid phases, attaining two different con- 
stant values in either phase, with a rapid transition region 
in the vicinity of the solidification front. The level set of 
(j>(r) — is identified with the solidification front, and the 
dynamics of (f> are carefully designed so that the level set 
dynamics follows that of the evolving solidification front 

The phase-field model finesses the problem of front 
tracking, but it is still prohibitively expensive for large 
systems. The cost is driven by the combined require- 
ments of fine resolution near the interface and a domain 
size set by the diffusion length and time for the system 
to evolve to steady state. The grid spacing must be 
small enough that the phase-field model converges to the 
Stefan problem, often referred to as the sharp interface 
limit. Recent work by Karma and Rappel p3| , involv- 
ing an improved representation for the temperature field 
within the interface, has extended this limit to the order 
of the capillary length, typically 10 _8 m. While improved 
asymptotics does help, the large system size needed to 
contain large diffusion-limited structures which form at 
low undercoolings, and the extended time scales have re- 
mained out of reach. Unfortunately, it is precisely these 
conditions that prevail in the most successful experiments 
done to date Q. 

Our contribution in this article is to show how Karma 
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and Rappel's phase field model can be implemented in a 
computationally efficient manner, thus removing a signif- 
icant obstacle to the numerical solution of large-scale so- 
lidification problems. We exploit the fact that the phase 
and temperature fields are both essentially constant over 
most of space, with significant variation only near the 
solidification front itself. This suggests that an adaptive 
mesh method which concentrates grid points in the inter- 
face region will be efficient. However, the implementation 
is non-trivial, and we have found it effective to use finite 
element methods, as described below. 

We first illustrate that our method allows faster and 
more efficient numerical integration of phase-field mod- 
els, especially in large systems and integration for long 
times. We then examine the physics of dendritic growth 
using a phase-field model solved by this method. In par- 
ticular, we present the convergence properties and veloc- 
ity selection of dendrites. We also measure the effective 
anisotropy introduced by our adapting grids, and com- 
pare it to that obtained using uniform grid methods. 

We model solidification using the phase-field model 
used by Karma and Rappel j3^|. We rescale tempera- 
ture T by U — cp(T — Tm)/L, where cp is the specific 
heat at constant pressure, L is the latent heat of fusion 
and Tm is the melting temperature. The order parame- 
ter is defined by <fi, with (ft = 1 in the solid, <f> — —1 in 
the liquid. The interface is defined by <f> — 0. We rescale 
time by r G , a time characterizing atomic movement in 
the interface, and length by W , a length characterizing 
the narrow region between liquid and solid. The model 
is given by 
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where D = ar /W 2 and a is the thermal diffusivity, and 
where A is a parameter that controls the coupling of U 
and (f>. Anisotropy has been introduced in Eqs. (Q) by 
defining the width of the interface to be W(n) — W A(n) 
and the characteristic time by r(n) = T A 2 (n) p3fl , with 
A(n) e [0, 1], and 
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The vector n in Eq. (2) is the normal to the contours of A, 
and 6 /x and 4> iV represent partial derivatives with respect 
to x and y. The constant e parameterizes the deviation 
of W(n) from W . 

Karma and Rappel |33| derived asymptotic relation- 
ships between the parameters of Eqs. (|l]) which allow 
them to operate in the sharp interface limit, where 
U at the interface satisfies t/i„t = —d(n)K — /3(n)V n , 



where d(n) is the capillary length, k is the local cur- 
vature, (3(n) is the interface attachment kinetic coeffi- 
cient and V n the normal speed of the interface, all as- 
sumed in dimensionless form here. In terms of Eq. (|^), 
d{n) = d [A(n) + d 2 A(n)] , where d = Ei/X, where 
Ei = 0.8839 [f33| and 6 is the angle between n and the 
cc-axis. Karma and Rappel |33| also showed that W , t 
and A can be chosen so as to simulate arbitrary values of 
/3. In particular, choosing A — C\D = Ciar /W 2 , with 
C\ = 1.5957, makes [3 = 0, a limit which is appropriate 
for SCN @. 

We compute four-fold symmetric dendrites in a 
quarter-infinite space, initiated by a small quarter disk 
of radius R centered at the origin. The order param- 
eter is initially set to its equilibrium value 4> (x) = 
— tanh((|s| — R )/V2) along the interface. The temper- 
ature is initialized to be everywhere equal to its far-field 
undercooling C/(|x|) = c^(Too — T m )/L = —A. 

We simulate Eqs. (Q) on an adaptive grid of lin- 
ear isoparametric quadrilateral and triangular finite el- 
ements, formulated using Galcrkin's method. Element 
data are arranged on a two dimensional clement-quadtree 
data structure |Q, making our code scalable. The grid 
is locally refined to have a higher density of elements in 
the vicinity of the interface of the <^-field, as well as in an 
extended region in the liquid which contains the [/-field. 
The criterion for refinement is based on changes in fluxes 
of both the 6 and U fields. Typically, the grid is adapted 
every 100 time steps which permits the 6 and U fields 
to remain within the refined range between regridding 
updates. We allow a difference of at most one level of 
refinement between neighboring quadrilateral elements. 
In such a case the quadrilateral element of lower level of 
refinement has an extra side node. The extra nodes are 
resolved with triangular elements. Details of our algo- 
rithm will be presented in an upcoming publication. 

When using an adaptive grid procedure, the concept of 
a grid spacing is replaced by that of a minimum grid spac- 
ing Ax m i n , representing the finest level of spatial resolu- 
tion. We found that best convergence is obtained when 
the algorithm layers the grid so the highest density of 
elements appears around the 6 interface, whose width is 
of order 1. The [/-field ahead of the A field is of order 
D/V n and has smaller gradients than 6, and so it is en- 
compassed by a uniform mesh of grid spacing 2Ax min 
and 4Aa; m ; n . We found that the convergence of our solu- 
tions is relatively insensitive to Ax m [ n . For a test case of 
dendrites grown at A = 0.55, D = 2, e = 0.05, and inte- 
gration time step dt — 0.016, our solutions for the steady 
state velocity converge to that given by solvability theory 
to within a few percent for 0.3 < Ax mnl < 1.6. 

Fig. shows a typical dendrite 10 5 times steps into 
its evolution computed using our adaptive grid method. 
For this case, A = 0.7, D = 2, e = 0.05, and the time 
step dt = 0.016. The system size is 800 x 800 with 
Ax m in = 0.8. Side branching is evident, and arises due 
to numerical noise. About half the computational do- 
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main is shown. This calculation took approximately 10 
cpu-hours on a Sun UltraSPARC 1200E workstation. 




FIG. 1. A dendrite grown using the adaptive-grid method 
for A = 0.55, D = 2, e = 0.05, and dt = 0.016. Clockwise, 
beginning at the upper right the figures show contours of the 
(7-field, the contour <j> — 0, contours of the (^-field and the 
current mesh. 



We examined the cpu-scalability of our adaptive grid 
algorithm with system size by growing dendrites in sys- 
tems of linear dimension L b and measuring the cpu time 
R1 for the dendrite to traverse the entire system. For 
A = 0.55, Ax m i n — 0.4 and the other parameters the 
same as those used in computing Fig.Jll, the relationship 
between and Lb is shown in Fig. gwhere we see that 
R c l ~ L 2 B . The number of calculations performed, per 
unit time, is proportional to the number of elements in 
our grid, which is set by the arclength of the dendritic 
interface being simulated 
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FIG. 2. CPU time vs. the system size, illustrating the 
quadratic dependence of computing time for a dendrite to 
move through the system on linear dimension Lb- 



multiplied by the diffusion length D/V n . For a 
parabolic shape the arclength is approximately Lb- 
Thus, since the dendrite moves at a constant velocity 
V n , 
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where is a constant that depends on the implemen- 
tation. The cpu time i?" needed to compute a full den- 
dritic microstructure on a uniform grid scales as R™ = 
[Rg/(y n Ax^ n )\L' B . For large system sizes, our method 
will be faster than uniform grid methods by a factor Lb- 
We tested the velocity selection mechanism of the 
phase-field model solved by our adaptive-grid algorithm 
for various undercoolings. In all cases we found very 
good agreement with the results of solvability theory. 
Fig. U shows the rescaled tip velocity V n d /D vs. time 
for three different undercoolings, A = 0.45, 0.55, 0.65, 
with corresponding dimensionless thermal diffusivities 
D = 3,2,1, and dimensionless capillary length d Q = 
0.185,0.277,0.554 respectively. In all cases the steady 
state tip velocity of our dendrite reproduces solvability 
theory - shown by the horizontal lines - within a few 
percent. 
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FIG. 3. The time evolution of the dimensionless tip velocity 

for undercooling A = 0.45, 0.55, 0.65. The horizontal lines are 

the results of solvability theory. 



We also tested the effective anisotropy of our dynami- 
cally adapting lattice following the method of Karma and 
Rappel (3^]. We couple <f> to an initially constant back- 
ground temperature = d /R . This maintains the 
interface in equilibrium for the case of isotropic surface 
energy. When growing in the presence of anisotropy, Ab 
is adjusted dynamically so as to maintain the velocity of 
the interface, measured along the x-axis, to zero. Even- 
tually our crystal neither shrinks nor grows, and its final 
equilibrium shape is fitted to the form 

R(9) = R o (l + e c «cos0) (4) 

where R is the equilibrium radial co-ordinate the crystal 
and 9 is the polar angle measured from its center. The ef- 
fective anisotropy represents the modification to e due to 
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the grid. Fig. ^ illustrates a crystal grown to equilibrium 
using an input anisotropy e = 0.04. Using to Eq. (^), we 
find e e ff = 0.041. This accuracy is typical. 
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FIG. 4. The equilibrium shape of the interface, for an input 
anisotropy e = 0.04. The effective anisotropy e = 0.041. 



This letter presents a new adaptive grid algorithm that 
is used to study solidification microstructures using adap- 
tive refinement of a finite element grid. Our method is 
used to solve the phase-field model given by Eqs. (Q). 
Our main result is that our solution time scales linearly 
with system size, rather than quadratically as one would 
expect in a uniform mesh. This allows us to solve the 
phase-field model in much larger systems and for longer 
simulation times. We showed that the convergence of our 
solutions remains accurate over a large range of Aa; m . 
Furthermore, dendritic tip speeds were found to be re- 
produced within a few percent of the theoretical values 
predicted by solvability theory. The effective anisotropy 
induced by our method was found to be a few percent of 
the input anisotropy. The speed increase of our method 
allows us to investigate dendritic microstructures at un- 
dercoolings somewhat lower than A = 0.1 in 2D. These 
results, as well as work in 3D will appear in an upcoming 
publication. 
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